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Abstract 

We develop a variational optimization method for crystal analysis in atomic 
resolution images, which uses information from a 2D synchrosqueezed trans¬ 
form (SST) as input. The synchrosqueezed transform is applied to extract 
initial information from atomic crystal images: crystal defects, rotations and 
the gradient of clastic deformation. The deformation gradient estimate is 
then improved outside the identified defect region via a variational approach, 
to obtain more robust results agreeing better with the physical constraints. 
The variational model is optimized by a nonlinear projected conjugate gradi¬ 
ent method. Both examples of images from computer simulations and imag¬ 
ing experiments are analyzed, with results demonstrating the effectiveness of 
the proposed method. 

Keywords: Atomic crystal images, crystal defect, clastic deformation, 2D 
synchrosqueezed transforms, optimization 


1. Introduction 

Defects, like dislocations, grain boundaries, vacancies, etc., play a funda¬ 
mental role in polycrystalline materials. They greatly change the material 
behavior from a perfect crystal and affect the macroscopic properties of the 
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materials. Analysis of images arising from atomistic simulations or imaging 
of polycrystalline materials hence becomes very important to characterize 
and help understand the defects and their effects on crystalline materials. 
While the defect analysis is traditionally done by visual inspection, the large 
amount of data made available due to advances in imaging and simulation 
techniques creates a need of efficient computer-assisted or automated analy¬ 
sis. 

Crystal deformation at the atomic scale is another important quantity 
that characterizes poly crystalline materials. When the deformation, denoted 
by 0, is well-defined, the tensor field F = V0 describes the local crystal 
strain; the polar decomposition of F at each point gives grain rotations; the 
curl of F provides information about defects and the well-known Burgers 
vector that represents the magnitude and direction of the lattice distortion 
resulting from a dislocation. Since it is almost impossible to estimate the 
deformation manually, the development of computer-aided analysis becomes 
important. 

For crystalline materials, defects are physical domains of the materials 
such that it is not possible to identify a smooth crystal deformation 0 = 0 -1 
that maps the atomic configuration back to a perfect lattice. In other words, 
the deformation gradient G = V0 = F is irregular and has nonzero curl at 
the defect location. In the opposite case, when a smooth deformation map 
does exist, the affine transform given by the gradient of the map, G = V0, 
transforms the image locally to an undistorted lattice of atoms. Therefore, 
for a defect-free region of the material, G is a gradient field and thus is 
curl-free: curlG = 0. Crystal image analysis hence requires the detection 
of the defect regions and preferably also the estimation of the local elastic 
deformation G away from the defects. 

In recent years several variational image processing methods for crystal 
analysis have emerged. The most basic versions just segment the crystal 
image into several crystal grains and identify their orientation, using clever 
formulations as convex optimization in 2D mm or efficient, sophisticated 
optimization algorithms for 3D images [9]. Berkels et al. additionally deter¬ 
mine a full deformation field 0 |2j via nonlinear optimization, while defects 
and local crystal distortion G are identified in [8], 

If accurate atom positions can easily be extracted from the image (or if 
the data stems from atomistic simulations) such that the input data consists 
of a discrete list of atom positions, then the local lattice orientation and 
deformation as well as defects can be obtained efficiently by identifying the 
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nearest neighbors of each atom mmm- 

Another efficient approach is the crystal image analysis via 2D syn- 
chrosqueezed transforms (SSTs) in [16]. The SST was originally proposed 
in [5], rigorously analyzed in [3J HI], [15], and extended to 2D in (TTl [19]. It is 
proved that the 2D SST can accurately estimate the local wave vector of a 
nonlinear wave-like image. Inspired by the fact that a deformed atomic crys¬ 
tal grain can be considered as a superposition of several nonlinear wave-like 
components, [16] proposes an efficient crystal image analysis method based 
on a 2D band-limited SST. In particular, tracking the irregularity of the syn- 
chrosqueezed energy can identify defects, and the deformation gradient G can 
be obtained by a linear system generated with local wave vector estimations. 
A recent paper [18] on the robustness of SSTs supports the application of 
this analysis method to noisy crystal images. The idea of exploiting the local 
periodic structure in the Fourier domain to extract a deformation gradient 
has also been considered in pTO] in the crystal imaging literature. 

Our work here is a variational approach based on the information obtained 
from a band-limited 2D synchrosqueezed wave packet transform following 
[16]. In this manuscript, we use this information as input to a variational 
optimization in order to improve the robustness of the analysis and, impor¬ 
tantly, to make the results better agree with the physical nature of defects. 

2. Variational method to retrieve deformation gradients 

Let us fix a perfect, unstrained crystal with a fixed orientation as our 
reference lattice. Let C R 2 be the domain of the image. Our objective is 
to fold at each igU the local strain or deformation gradient G(x ) G R 2x2 of 
a (locally defined) deformation </> which deforms a defect-free neighborhood 
of x into the reference lattice. (Note that since the image corresponds to a 
deformed crystal state, x should be understood as an Eulerian coordinate.) 
In other words, we seek the affine map defined by G(x) which maps the local 
atom arrangement to the reference lattice. This is however impossible around 
defects. In particular, while curl G — 0 in the elastic region, in the defect 
region, curl G is not zero, and the integral of curl G over a neighborhood of 
the defect such as a dislocation should match the defect’s Burgers vector [[9]. 

Assume the defect region is given by C D and the curl of G is given 
by b, consistent with the Burgers vectors and with curl G = 0 on \ We 
expect the displacement held to minimize the elastic energy of the system 
outside the defect region, since the system under imaging is in a quasistatic 
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state. Given G 0 a rough guess of the deformation gradient, this motivates 
the energy minimization 

min [ \G-G 0 \ 2 + W(G)dx 
G Jn\n d (1) 

s. t. curl G = b , 

where |-| denotes the Frobenius norm of a matrix, |A| = (tr(A T A)) 1//2 , and 
W is the elastic stored energy density. 

Since our reference lattice represents the undeformed equilibrium state of 
the crystal and the atom configuration in the image is produced by the (local) 
deformation 0 _1 , the stored elastic energy can be expressed in the standard 
Lagrangian form as the integral over the reference domain 0(fi \ Qj) of an 
elastic energy density w that depends on V(^)^ 1 ) = G o 0 -1 , 

[ w(G~ l of^t/)) dy. 

J^n\n d ) 

Here, w satishes the standard conditions coming from hrst principles, i. e. w 
is frame indifferent, w(A) = 0 for A e SO( 2), w(A) > 0 else, and w(A) = oo 
if det A < 0. After a change of variables the elastic energy turns into 



W(G) dx 


for W (G) = w{G~ x ) det G, where it is easy to see that W has the same above 
properties as w. For w (or equivalently W) one can use a material-specific, 
possibly anisotropic energy density. To be specific, since our numerical ex¬ 
amples are all concerned with a triangular lattice exhibiting isotropic elastic 
behavior, we here simply restrict ourselves to the following neo-Hookean-type 
elastic energy density 

W(G) = |(|G| 2 - 2) + (| + \)(det G - l) 2 (2) 

— /r(det G — 1). 


Note that in (JTJ) , the fidelity and elastic energy terms are both evaluated 
outside the defect region. Within the defect region, since it is not possible to 
map the local configuration of atoms back to the reference state, the estimate 
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G 0 is not trustworthy. It is also well known that the elastic energy blows up 
logarithmically approaching the dislocation core, and hence it only makes 
sense to penalize the elastic energy away from the defects. 

The questions are then how to estimate G 0 and to determine the defect 
region as well as the vector held b consistent with the defects’ Burgers 
vectors. In this work, the desired information is all obtained using syn- 
chrosqueezed wave packet transforms. We first review the synchrosqueezed 
wave packet transforms and the estimate of Go in §2.1[ The defect region 
and Burgers vectors are estimated using the synchrosqueezed wave packet 
transforms as explained in §2.2 The variational problem (JTj) is then solved 
by a constrained minimization algorithm as described in 


2.1. Synchrosqueezed wave packet transforms (SSTs) 

Before explaining the variational optimization, we briefly introduce the 
crystal image model established in [TB] and the synchrosqueezed wave packet 
transforms in [171 IP] for the purpose of a self-contained description. Consider 
an image of a polycrystalline material with atomic resolution. Denote the 
perfect reference lattice as 


C = {av i + bv 2 : a,b integers} , 

where v \, G M 2 represent two fixed lattice vectors. Let s(a:) be the image 
corresponding to a single perfect unit cell, extended periodically in x with 
respect to the reference crystal lattice. We denote by D the domain occupied 
by the whole image and by 12*,, k — 1,..., M, the grains the system consists 
of. Now we model a polycrystal image / : D —y M as 


f(x) = a k (x)s(N(p k (x)) + c k (x) if x e fl k , (3) 

where N is the reciprocal lattice parameter (or rather the relative reciprocal 
lattice parameter as we will normalize the dimension of the image). The 
<f k : Vt k —y M 2 is chosen to map the atoms of grain Q k back to the configuration 
of a perfect crystal, in other words, it can be thought of as the inverse of 
the elastic displacement held. The local inverse deformation gradient is then 
given by G = in each D k . For generality, ^ also includes a smooth 
amplitude envelope a k (x) and a smooth trend function c k (x) to take into 
account possible variation of intensity, illumination, etc. during the imaging 
process. Using the 2D Fourier series s of s and the indicator functions Xn k , 
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we can rewrite (j3]) as 


M 

f ( x ) = 

k =1 


,£e£* 


+ Cfc(x) 


( 4 ) 


where £* is the reciprocal lattice of C (recall that s is periodic with respect 
to the lattice £). In each grain 12*,, the image is a superposition of wave-like 
components ak{x)'s{^)e lN ^^ k ^ with local wave vectors iVV(£•</>&(a;)) and lo¬ 
cal amplitude a/ c (:r)|s'(£)|. Our goal here is to apply the band-limited 2D SST 
developed in [16j to estimate the defect region and also Go = J2i=iXn k 'V4 , k 
in the interior of each grain Q*. as required in the variational formulation. 

The starting point of 2D SST is a wave packet w a o x , which is, roughly 
speaking, obtained by translating, rotating, and rescaling or modulating a 
mother wave packet w : M 2 —» C according to the parameters x G M 2 , 
0 G [0,27r), and a G M [Hi 331 HE]- Let Wf(a,9,x) = j n w^.{y)f(y) d y be 
the wave packet transform of / at scale a and angle 6 in the frequency domain, 
and at spatial location x. The wave packet transform is a generalization of 
curvelet and wavelet transforms with better flexibility in frequency scaling 
and consequently is better suited to analyze crystal images with complex 
geometry. As a convolution with smooth wave packets, Wf is well-defined 
and smooth even under very low regularity requirements for /, e. g. / G 
L°°(M 2 ). The spectrum of the windowed Fourier transform of a given crystal 
image spreads out in the phase space, as illustrated in Figure [ljb) . The 
2D synchrosqueezed wave packet transform (SST) aims at sharpening this 
phase space representation. In the SST, for each ( a,9,x ), we define the 
corresponding local wave vector 


Vf(a, 9 , x) 


2niWf(a, 9, x) 


Here, V x Wf denotes the gradient of Wf with respect to its third argument 
x. The synchrosqueezed (SS) energy distribution of / is then constructed as 


Tf(v,x)= / / \Wf(a,9,x)\ 2 5(vf(a,9,x) — v)adad9, (5) 

Jo Jo 

where 6 denotes the Dirac measure. By the stationary phase method, in a 
small neighborhood of a point x such that Wf(a,9,x) is nonzero for some 
fixed (a, 9), Wf(a,9,x) is essentially a plane wave in x with a local wave 
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vector approximating a certain WV(£ • <pk{x))- Inspired by this intuition, it 
is proved in jl6] that in the interior of a grain, the local wave vector estimation 
Vf(a,9,x ) approximates local wave vectors WV(£ • <j>k(x)) accurately, if the 
deformation £>*, and the amplitude function are smooth. The SST squeezes 
the wave packet spectrum \Wf(a,9,x)\ 2 according to Vf(a, 9,x) to obtain a 
sharpened and concentrated representation of the image in the phase space. 
Hence, in the interior of a grain, the SS energy distribution Tf has a support 
concentrating around local wave vectors IVV(£ • cftk(x)), £ G £*, and is given 
approximately by (see e.g., Figure [lje) in polar coordinates) 

T f (v,x) « ^a fe (x) 2 |s(O| 2 5(u«7VV(e-0 fc (a;))), (6) 

£e£* 

understood in the distributional sense. Therefore, by locating the energy 
peaks of Tf, we can obtain estimates of local wave vectors 7VV(£ • (j>k{x)) and 
also their associated spectral energy. In practice, we will aim at high energy 
peaks corresponding to £ close to the origin in the reciprocal lattice. For 
simplicity and concreteness, we will in this article focus on the case when 
the lattice is hexagonal. We then have six such reciprocal lattice vectors 
£, which can be further reduced to three due to the symmetry £ gg —£. 
We will henceforth denote these as £ n , n = 1,2,3, and denote by u® st (x) 
the estimate of lVV(£ n • (j>k{x)) = N(V(f>k(x))£ n . The inverse deformation 
gradient G 0 (x) = V(j>k(x) is determined by a least squares fitting of the 
deformed reference reciprocal lattice vectors lV£ n to the u® st : 

3 

G 0 (a;) = argmin V ||< st (a:) - NG£ n \\ 2 . 

G 

In practice, for each physical point x we represent Tf(-,x) in polar coordinates 
(r, $) G [0, oo) x [0,7r) (where $ G [vr, 27 t) is ignored due to symmetry). To 
identify the peak locations {u® st }, we choose the grid point with highest 
amplitude in each 60 degree sector of d. 

As the representation © of Tf is no longer valid around the crystal 
defects, we may characterize the defect region by using an indicator function 
for the deviation of Tf from the representation ([6]). To this end, for each 
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n G {1, 2, 3} (corresponding to one of the sectors), we define 


w n (x) 



est\ 
n ) 


Tf(v, x ) dn 


' argi;G[(n—l)7r/3,n7r/3) 


Tf(v, b ) dv 


where Bfiv^f') denotes a small ball around the estimated local wave vector 
v® st . Hence, rnass(x) := 'Yh n w n{x) will be close to 3 in the interior of a grain 
due to © , while its value will be much smaller than 3 near the defects. This 
is illustrated in Figure [ljc) , where we show rnass(x) for the crystal image in 
Figure [l](a) . The estimate of defect regions can be obtained by a thresholding 
mass(x) at some value rj G (0, 3) according to 


fii = {i G H : mass(x) < g} , 

an illustration of which is shown in Figure [l](d) . Figure [2] shows the estimate 
of G 0 by the SST. It is observed that the crystal is slightly deformed in the 
grain interior and heavily deformed at the defect region. 

To better interpret the inverse deformation gradient Go, we compute its 
polar decomposition G 0 (x) = U 0 (x)P 0 (x) for each point iGfl, where U 0 (x) is 
a rotation matrix and Pq(x) is a positive-semidehnite symmetric matrix. The 
rotation angle of Ufix) describes the crystal orientation at x: det(Go(x)) — 1 
indicates the volume distortion of Go(x)] the quantity |Ai(rr) — A 2 (x)|, where 
Ai(x) and A 2 (a;) are the eigenvalues of -Po^), characterizes the difference in 
the principal stretches of G 0 (ic) as a measure of shear strength. The bottom 
panel of Figure [2] shows these quantities corresponding to the estimate of 
Go in the top panel. In the later numerical examples, instead of G itself we 
will always present the crystal orientation, the volume distortion, and the 
difference in principal stretches. 


2.2. Defect region, topological consistency, and Burgers vector identification 
Based on the estimated Go, we can determine the defect region and 
also the Burgers vectors or equivalently the curl held b used in the constraint 
of our variational formulation ([!]). We give the details in this section. 

2.2.1. Burgers vectors and curlG 

As explained previously, away from defects, G can be interpreted as the 
gradient V0 of a (fictitious) deformation 0 deforming the configuration of the 





Figure 1: (a) An example of a crystal image, (b) Windowed Fourier transform at a local 
patch indicated by a rectangle, (c) The defect indicator mass(a;) provided by the SST. 
(d) Identified defect region for a threshold of rj = 2. (e) The SS energy distribution in 
polar coordinates at a point outside the defect region, (f) SS energy distribution in polar 
coordinates at a point in the defect region. 


given image into a perfect reference crystal of a fixed orientation. Now the 
fact that gradient fields are always curl-free can be exploited as a constraint 

curl G = 0 on Q \ Q f i (7) 

in the variational method. In the defect region Qj, however, the interpreta¬ 
tion of G as a deformation gradient breaks down since there is no smooth 
deformation of the crystal that can undo the lattice defect. In fact, denoting 
the connected components of f Id by O^,..., Q l d , it is relatively simple to see 
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Volume Distortion 


Figure 2: Top panel: Estimated inverse deformation gradient Go € R 2x2 of the atomic 
crystal image in Figure [I] (a). Bottom panel: The crystal orientation, the difference in 
principal stretches, and the volume distortion of Go- The grey mask in these figures is the 
defect region identified in Figure [l] (d). 

that the integral 

B: t = j curl G dx (8) 


10 























Figure 3: A closed curve 7 around a defect in the crystal image breaks up when transported 
to the perfect crystal, the gap being the Burgers vector. 


is the Burgers vector associated with the defect in Q d . Indeed, consider 
any closed injective curve 7 : [ 0 , 1 ] —* M 2 that encloses Q d counterclockwise, 
but does not intersect Figure [3j left, shows a specific example of such a 
curve, which is chosen such that it connects a sequence of neighboring atoms. 
Denote the curve interior by fl. Since dfl lies in the defect-free region, for 
every x € <9f2 there is a deformation (j) (the gradient of which is given by 
G) that deforms a neighborhood of x into the perfect reference configuration 
of the crystal. Starting at 7 ( 0 ) we can thus iteratively map little segments 
of the curve 7 into the reference crystal via 7 o 0 -1 , resulting in a curve 
7 : [0,1] -A M 2 on the reference configuration. Note that 7 is in general no 
longer closed or injective. By Stokes’ Theorem, using the tangent vector n L 
to dfl, we have 


curl G dx = / curl G dx = 


Gn ± da 


Ian 


= G('Y(t))i(t)dt= 7 (f) dt = 7 ( 1 ) — 7 ( 0 ). 


The vector on the right-hand side is the defect’s Burgers vector and is obvi¬ 
ously independent of the originally chosen curve 7 . If Q l d contains an isolated 
dislocation surrounded by a regular lattice, B { just represents the Burgers 
vector of that dislocation. If Q d contains multiple dislocations or even a sec¬ 
tion of a high angle grain boundary, Bi represents the accumulated Burgers 
vector of all defects in Vt l d (note that a high angle grain boundary may be 
thought of as a string of dislocations with distance smaller than the lattice 
spacing so that the single dislocations are not clearly spatially separated; 
Figure | 6 ]( c) gives an impression of the size of the numerically found Q l d . 

As in the case of \ f^, where we know curl G = 0 , we also have a priori 
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information on curl G in fl l d . In particular we know that Bi is a Burgers 
vector and thus must lie in the discrete set of Bravais lattice vectors of the 
perfect reference lattice C. We thus identify Bi by projecting the (potentially 
noisy) estimate f Qi curl G 0 dx onto C and then impose ([8]) as a constraint on 
G. In fact, instead of prescribing the accumulated curl in fl d via (|8]) we may 
just as well prescribe 

curl G = bi on fl d (9) 

for a function bi : fl l d —> R 2 with f Qi bi dx = Bi (in mathematically more 
precise terms, bi may be a distribution). This is possible since we are only 
interested in the held G outside of fld and since any held G : fl —> R 2x2 
satisfying ([7]) and ([8]) can be modified on fld to a held satisfying ([9]). The 
function bi is here simply chosen as bi = diag(a, j3) curl G 0 with a, /3 G R 
such that f Qi bi dx = Bi (he., an overall scaling of curl G' 0 ). Summarizing, in 
our variational method to extract G from the noisy data Go we will prescribe 
the constraint 


curl G = b for b 


0 on Q \ fl d , 
bi on fl l d . 


( 10 ) 


2.2.2. Refined defect regions 

On the one hand, the threshold r] to identify fid should be chosen very 
low to yield thin and localized defect regions (e. g. such that defect regions 
fl d around single dislocations stay separated from each other), on the other 
hand, the thinner the identified defect region f l d the worse will the estimate 
of the Burgers vectors Bi be. A compromise is to first use thick defect 
regions fl l d in order to estimate the Burgers vectors Bi and then to impose 
the constraint (10) with a much finer estimate of the fl d . However, it may 
happen that a thick patch fl l d contains multiple thin connected components 
flf,..., flf. In that case the B tj are defined as the closest projections of 
f ij curl G 0 dx onto C under the constraint B ix + ... + B ik = Bi, where Bi is 

a 


the accumulated Burgers vector of the patch fl d . In order to obtain a very 
thin and localized fld we simply identify the ridge of 3 — mass(6) inside the 
thick fld and then dilate this ridge by a few pixels. The ridge computation 
turns out to be sufficiently robust even in the presence of noise, as can be 
seen from the similarity between the detected defect regions with and without 
noise in Figures [6] and [7j 
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Figure 4: Along a closed path 7 traversing a sequence of crystal grains, the deformation 
gradient G changes continuously from I to R 7 ^ I. The gray shade indicates the local 
crystal orientation from the identity I (white) to R (dark gray). Dots represent point 
dislocations; lines indicate high angle grain boundaries. Along the path 7 all grains are 
connected by low angle grain boundaries. 


2.2.3. Introduction of jump sets for topological consistency 

A Bravais lattice does typically not only exhibit translational, but also 
rotational symmetry. The so-called point group P C SO( 2) comprises all 
those rotations which leave the reference lattice invariant. This leads to an 
ambiguity in the deformation gradient G: if G correctly describes the local 
configuration of the crystal, then RG for any R e P does so as well. Even 


though the constraint (10) has the effect that the matrix field G will be 


locally consistent (in the sense that for any y in a neighborhood of x G Cl we 
have \G(y) — G(x)\ = min Rg p \RG(y) — G(x)\), global consistency is often not 
guaranteed. Indeed, Figure [4] shows a situation in which along a closed path 
7 C 11, G changes continuously from the identity / to an element R 7 ^ / of 
the point group. Since R describes the same local crystal configuration as /, 
the curl where G jumps from R to I is spurious. As in |Sj, this inconsistency 
can be remedied by introducing a cut set S across which the tangential 
component of G is allowed to jump by a point group element, 


G n L = RG + n ± 


for some R : S -A P, where G~ and G + denote the value of G on either side 
of S and n 1 - : S —> M 2 denotes the tangent to S. 

Note that we have large freedom to choose the cut set S. For instance 
we could take S to be composed of smooth lines connecting the different 
connected components Q l d of fid- Those lines can easily be chosen in such a 
way that they do not intersect and that the connected components, say IT, 
of f] \ (fid U S ) are simply connected. Within each IT, the crystal is defect- 
free, thus there is a deformation 0 * mapping IT onto the reference crystal. 
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The matrix field G will be the gradient of 0* inside fF, which is automati¬ 
cally consistent all over fib If G jumps across S between two neighboring 
components fF and fF, then for all x € S '*- 7 := <9fF D <9fF we must have 
G t (x)n ± (x ) = RG^ (x)n- 1 (x) with a constant R e P, since S * J lies inside a 
defect-free crystal region (G* and G - 7 denote the value of G in fF and fF, 
respectively, while n L is the tangent vector to S ). Therefore, not only is S 
smooth, but even R constant on each curve segment of S. Furthermore, up 
to multiplications with point group elements, the optimal G does not de¬ 
pend on the particular position and geometry of S. Indeed, if S '*- 7 is shifted, 
producing new IF, fF, S'*- 7 , then the new optimal field G is given by 

{ G(x) if x G IF fl fF or x G fF fl fF, 

RG{x ) if x e fF n fF, 

R~ 1 G(x) if x e fF fl Q*. 

For simplicity, the cut set S' in our computation is algorithmically chosen 
after the problem discretization. First we sweep over all image pixels and 
reassign the pixel values via a fast marching type method. Let Q denote the 
set of already treated pixels (consisting initially of only one pixel). Amongst 
the not yet treated neighbors of all pixels in Q , we pick the pixel x with 
the smallest value of miiiRgp |/2Go(x) — Go(x n )|, where x n G Q denotes the 
neighboring pixel. We replace its value Go(x) with RGq(x) for the minimizer 
R £ P and set Q = Q U {x}. After all pixels have been swept through in 
this way, we set S to be the union of interfaces between neighboring pixels 
Xi,X 2 that have min Rg p |i?G 0 (xi) — G 0 (x 2 )| < |G 0 (xi) — G 0 (x 2 )|. Note that 
we also take care that Qd D S = 0 so that the estimation of Burgers vectors 
is not impaired. 

2.2.4 ■ Variational principle with topological jump set 

While we will use G on the whole domain in the numerical algorithm, 
as the energy (|TJ) depends only on G on Q \ Q d , it is more convenient to 
reformulate the problem with a bounded Lipschitz domain f2 0 := f2\ (f^US 1 ) 
and consider the Hilbert space 

if(curl,ft 0 ) = [ue (L 2 (n 0 )) 2 

with norm given by 

IMItf(curlA)) = ll' U ll(L 2 (O 0 )) 2 + H c Urlt(|| L 2 ( n o) . (12) 


curlu G L 2 (f2 0 )|, (11) 
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Note that dfl 0 — S U dfld = S U jj. dtt d , where we recall that fl d are the 
connected components of Cld. We define the tangential trace operators 7 t _ 
and 7 + on the two sides of S for u G (C' oo (r2 0 )) 2 as 

It ( u ) = n± ■ u \s- and 7 t + (tt) = n L • w| 5 +, (13) 

where n = (rti, n-fi is the unit normal of S and vfi = (— ri 2 , nfi). It is standard 
(see e. g., i) using Green formula that 7 f can be extended to bounded 
linear operators from //(curl, Q 0 ) to the fractional Hilbert space H~ l / 2 (S). 
Similarly, we may also define the trace operator on dfl d . 

For the (inverse) deformation gradient G, it is then natural to consider 
the space //(curl, Q 0 ) 2 as the rows of G: G 1 ' and G 2 lie in //’(curl, h2 0 )- 
Thus, we generalize the tangential trace operators 7 ^ and 7 1 also to G, such 
that they act on each row of G and map //(curl, fi 0 ) 2 to {H~ l / 2 (S)) 2 and 
(i/ _ 1 // 2 (<9fIrf)) 2 , respectively. In particular, as we are only concerned with 
G on f2 0 , the constraint curlG = b on is equivalent to f gQi jfiG) da = 
f n bi dx. 

a _ 

With this, we may formulate (I1J) with the topological jump set more 
precisely as 


mm 


\G-G 0 \ 2 + W(G) dx 


Geif(curi,n 0 ) 2 Ja\n d 

s. t. curl G = 0 on fl 0 


7t(G) da = / bi dx 


'dn\ 


'n* 


1 ~{G) = R 1 t{G) on S, 


(14) 


for some fixed piecewise constant R : S —$■ P. Here the last equation 7 fi(G) = 
/? 7 t + (G) is understood in the space of (H^ l ^ 2 {S)) 2 (this is well defined as 
Rf G (// 1 / 2 (S ')) 2 for any / G (H l / 2 (S)) 2 due to the regularity of R). 


Theorem 1 . Assume that the topological jump set S and the boundary of fid 
are Lipschitz. Given Go G Al2x2(-h 2 (fA^))? R '■ S —>■ P piecewise constant, 
and b a finite measure on fid, the minimizer to the variational problem (14) 
exists. 
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Proof. We consider the admissible space of G: 


A — |g E H{ curl, f2 c 


curlG = 0 on O 0; 




lt{G) dc = / 6* da;, 

Jw d 

lt{G) = i*y+(G) on S k 


Note that this is a closed linear subspace of H( curl, f2o) 2 , and hence any 
bounded sequence is weakly compact in A. Given a minimizing se¬ 


quence of (14), as the energy is uniformly bounded, we have 


/ \G (n) - G 0 \ 2 dx < C 

for some constant C. Together with curl G-"-* = 0, we obtain that ||Gkd ||^( curl 

is uniformly bounded, and hence weakly convergent (in //(curl, f2o) 2 ) to 
G(°°) e *4.. The existence then follows from the lower semi-continuity of 
J|G^ — G 0 | 2 dx and / W(G) da; with respect to the weak (L 2 (f2 0 \ fl d )) 2x2 
topology (which is clearly weaker than the weak topology of //’(curl, f^o) 2 ). 

□ 


2.3. Constrained minimization algorithm 

In this section we describe the numerical algorithm to solve the minimiza¬ 
tion problem ({TJ) or rather the corresponding version with topological jump 
set S (14). It is more convenient to describe the algorithm for the following 
simplified formulation, which is equivalent after discretization. 


min [ |G-Go| 2 +W(G)da; (15) 

G:M2X2 J^ d 

s. t. curl G = b on \ S , G~ = RG + on S , 

where G~ and G + denote the value of G on either side of S. For simplicity 
we use the same notation for the continuous and discrete objects. 


2.3.1. Finite difference discretization 

The image domain is discretized via M ■ N Cartesian pixels, indexed 
by o;-;?/-position (m, n) £ = {1,..., M } x {1,..., N} (the pixel spacing is 

assumed to be one). For an index m we denote by m + = m + 1 and m~ = 
m — 1 the next larger or next smaller index, where for simplicity we assume 
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periodic boundary conditions and use cyclic indexing, i. e. (. M + ,n ) = (1 ,n), 
( m,N + ) = (m, 1), (l _ ,n) = ( M,n ), (m, 1 _ ) = (m,N) for all (■ m,n ) G fh 
The matrix fields are discretized accordingly, G, Go G (M 2x2 ) MxA b The jump 
set S follows the edges between the pixels and is represented as a collection 
of horizontal or vertical pixel pairs, S C {((m, n), ( k,l )) G x : (k,l) — 
(■ m + , n ) or (k, l ) = (m, n + )}. Furthermore, we define the function R : S —>• P 
such that R(m,n),(k,i) is the point group element with smallest distance to 
(Go)m,n(Go)'jj j- R is extended to b2 x Q by the identity in P. Derivatives 
in x- and ^/-direction are replaced by finite differences that respect the point 
group equivalence across S, 

- G«„, 

(^G)« „ = (fl (ra ,„), (ra ,„ + )G ra ,„ + )« - 0«,, 

where superscript ij denotes the (i,j)-matrix entry. Throughout, a super¬ 
script d denotes discrete differential operators. In particular, the discrete 
curl and Laplacian are defined as 

cmf* : (M 2 x 2 )MxJV ^ ^MxN ^ 

(curl d G) m>n = d d G'^ n - d d G'^ n , 

A d : (R 2 ) MxN -> (R 2 ) MxN , 

A d = — curl d (curl d )* , 


where G'jj, n denotes the i th column of the matrix G m , n and the superscript * 
denotes the adjoint operator, which in this particular case is given by 


2 \MxN 


(curb 1 )* : ^ ) 
((curl d )*I/) m;n = 




0 2x2 \MxN 


t / _ r>T t / 

v m,n lx (m,n-),(m,n) y m,n- 


tdT t/ _ T/" 

- rL (?ra-,n),(m,n) v m~,n v m,n 


2.3.2. Projected descent in constraint space 
The discrete optimization problem reads 


mincec* E[G\ 

for E[G\ = Yl(m,n)e n i\^m,n ~ ( G 0 )m,n\ 2 + W(G m! n)) 
and C b = {Ge (R^ x2 ) MxN : curl d G = b} . 
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The constraint space Cf, is an affine space and can be expressed as G + Cq for 
a G with cur]' / G = b. Hence, the energy can be minimized using a standard 
projected nonlinear conjugate gradient (NCG) descent in C &. In more detail, 
we employ a Fletcher-Reeves NCG method in which the derivative 8qE of 
E with respect to G is always orthogonally projected onto C 0 (i. e. onto its 
component parallel to Cb) so that the algorithm is performed within the sub¬ 
space C h . Due to accumulating numerical errors we also have to project the 
current estimate G back onto Cb from time to time. Denoting the projection 
onto Cb by proj c , the NCG algorithm is initialized with G = proj G Go- 
The projection proj c F is the solution to the constraint minimization 
ni in curl << c=b E m n I ^rn,n ~ F mjn \ 2 , which satisfies the optimality conditions 

b = curl' 1 F, 0 = F - G + (curl d )*A 

for a Lagrange multiplier A G (W 2 ) MxN . Applying cuiT z to the second equa¬ 
tion we obtain 

curl d G — b— - A d A. 

Note that ker(A d ) _L range(cuiT , j. Denoting by (—A d ) _1 : range(A d ) —y 
ker(A d ) ± the inverse of —A d , we obtain A = (—A d ) _1 (curl <i G — b) and thus 

proj^G = F = G — (curl cZ )*(—A d ) _1 (cuiT z G — b ). 

Once G is computed, the projection onto Cb can also be obtained as 

proj C( F = G + proj Co (F - G). 

2.3.3. Fast projection onto constraint space 

For the above projection we need to invert the discrete Laplacian opera¬ 
tor. Using periodic boundary conditions this would be very fast using FFT if 
there was no jump set S. However, with nonempty jump set S, our finite dif¬ 
ference operators do not turn into pointwise multiplications in Fourier space. 
In order to obtain a fast inversion we decompose A d = Aq + J, where Aq is 
the standard discrete Laplacian and J the linear operator accounting for the 
point group, 


A d ■ 

ZA 0 . 


d 2 \MxN 


-A 


t>2\MxN 


(Aq V)m,n — Vm-,n + Kn+,n + Vm,n- + Kn ,71+ 4V^77 , j7 t, ? 
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J : 

(R 2 ) MxN 

(M 2 ) Mxj 

= 


-I)V m 

+ 


I)v mr 

+ 


I)v m , 

+ 


- I)V„, 


Note that J is symmetric, and it is highly sparse and thus has a very small 
range L = range J and a large kernel ker J = ker J 1 = L l . If we decompose 
V E (R 2 ) MxN into 

V — Vl + Vi± G 1/ © L 1 , 

then we obtain 


-A d V = B -AqV l ± — (Aq + J)V L + B . (16) 

The solvability condition tells us that ker Aq is orthogonal to the right-hand 
side, so we can just as well project the right-hand side onto (ker Aq)- 1 by 
subtracting the mean, 

-A2Ux = proj (kerA rf,j. ((A? + J)V L + B) . 

Now choosing a basis L = (iq,..., vk) and ker Aq = (si, s 2 ), we write 

K 

1 /. ^ ^ AjUj, 

*=1 


Vtd — °i s i + a 2&2 + (— ■ Aq) 1 proj (ker ((Aq + J)Vl + B) 
— aisi + & 2 S 2 + (—Aq) 1 proj( ker 

K 

+ Xi (( _A o) _ 1 P r °j (ke r A^JVi ~ a) , 


where (—Aq) - 1 : (ker Aq)^ —y (ker Aq)- 1 denotes the inverse of — Aq. The 
degrees of freedom Ai,..., A k, cq, a 2 now have to satisfy the K + 2 equations 
Vj ■ V L ± — 0, j — 1,..., K , (due to Vj G L ) and s 3 ■ ((Aq + J)Vl + B) = 0, 
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j = 1,2, (the solvability condition for (fT6|) where the dot denotes the dot 
product. In detail, the equations are given by 

0 = diVj • Si + d 2 Vj ■ S 2 + Vj ■ ( —Ag)” 1 proj(ker 1 

+ E*=1 A * (((- A o) _1 P r °j(ker -JVi- Vj V^j , l (17) 

0 = Sj ■ B + J2i=l ^i s j ■ JVi . J 

To solve for A := (Ai,..., A k, «i, fl 2 ) T , we write the linear system as a matrix 
vector equation. Since we need to solve systems of the form — A d V = B 
several times for different values of B , we perform a QR decomposition, 

QRPX = rhs, 


where P is a permutation matrix such that the lower rows of R are zero 
(this is possible since the solution to — A d V = B is only uniquely specified 
up to a two-dimensional subspace due to dim ker A d = 2) and all other 
rows have nonzero diagonal elements. Note that this decomposition is the 
bottleneck of the algorithm with a complexity of 0(K 3 ). The degrees of 
freedom corresponding to the zero-rows can now be chosen freely (say, as 
zero), and the resulting 


V 


cqsi + a 2 s 2 + (-Aq) Voj (kerA d)xR 

K 

+ A *(- A o) _1 P r °j {ke r AgR J y i 
2—1 


r k 


cqsi + a 2 s 2 + (—Aq) 1 proj( ker A g)^- 


A* Jvi + B 

. £=1 


solves B = —A d V. Note that (—Aq) 1 can readily be computed via FFT and 
that the Vi can be chosen as shifted versions of v J G (M 2 ) MxAr , j = 1, 2, with 
v 3 = 0 except for the j th entry of v[^ being one. Thus, (—AQ) _1 proj( ker A d)±Uj 
can be evaluated efficiently as just a shift of (—Ag) _1 proj ( ' ker a^iF. 


3. Numerical results 

In this section, we present several numerical examples for images com¬ 
ing from both computer simulations and real experiments to illustrate the 
performance of our method. The corresponding code is open source and 
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available as SynCrystal at https://github.com/SynCrystal/SynCrystal, 
We first apply our method to analyze synthetic crystal images of phase held 
crystals (PFC) To show the robustness of our method, both noiseless 
and noisy examples of PFC images are presented. Moreover, we examine 
different kinds of crystal images from real experiments: (1) a TEM-image 
of GaN with a series of isolated defects that leads to a large angle grain 
boundary; (2) a photograph of a bubble raft with a thick transition region 
at grain boundaries; (3) a TEM-image of A1 with a noisy and irregular grain 
boundary. 

In all of these examples, we compare the initial estimated strain Go pro¬ 
vided by the SST-based analysis described in (2T and the improved results 
G from the variational method (15), where each time we display crystal ori¬ 
entations, difference in principal stretches, and volume distortion. For better 
visualization we mask out the identified defect regions (in which there is no 
meaningful notion of strain). The curl of Go (which in general violates the 
physical constraint of being zero outside fid and of being compatible with 
the defects’ Burgers vectors) as well as the average curl G per connected de¬ 
fect region (which is compatible with the defects’ Burgers vectors) are also 
shown. 

The main parameters for the SST-based analysis are two geometric scaling 
parameters s and t in the 2D SST (for details see [TH|) - Smaller scaling 
parameters result in better robustness of SSTs while larger scaling parameters 
give more accurate estimates in noiseless cases ra- Hence, we adopt t — s ~ 
1 in the examples with less noise and use t — s ~ 0.8 when the crystal image is 
noisy. As discussed in [18], for images with heavy noise, the synchrosqueezed 
transform can still provide reasonable initial guess via a highly redundant 
transform with more computational cost. The variational model parameters 
A and /j in (J2| are simply set to 1 in all of the following examples. 

During the synchrosqueezed transform we also downsample the original 
image in the frequency domain and reduce the mesh size for the variational 
optimization. As a result, the computational mesh of the optimization model 
is actually independent of the resolution of the crystal images. To achieve 
a short computation time, we typically choose the downsampling rate as 
large as possible such that defects can still be localized with a precision 
below one atom spacing. Figure[8] will show a comparison between different 
downsampling rates. 

The main computational cost arises from inverting the Laplace operator 
via the solution of 0 - The number of variables is proportional to the 
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Figure 5: A noise-free PFC image (left) and its noisy version (right) with a zoom-in 
detailing the marked part. 

number K of pixel pairs of the singularity set S, and the solution cost scales 
like A' 3 . Note that, since S has codimension 1, K behaves roughly like the 
number of image pixels along one direction (and not like the total number 
of pixels), so the computations are still feasible for quite large images. The 
runtime for the examples in Figure[9]is less than 10 seconds, the runtime for 
the 1024 x 1024 PFC image in Figure[5j left, is 20 seconds. The runtime for 
the example in Figure[l2] with heavy noise is about 1 minute due to extra 
effort to obtain robustness of the synchrosqueezed transform. 

3.1. Synthetic crystal images 

The PFC model is a well-established method to simulate elastic and plas¬ 
tic deformations, free surfaces, and multiple crystal orientations in nonequi¬ 
librium processes. Its simulation can provide synthetic crystal images with 
various topological defects. A noiseless PFC image is given in Figure [5] (left). 
It contains isolated defects and large as well as small angle grain boundaries, 
the latter showing up as a string of dislocations. 

Compared to the initial strain estimate G o, the optimized G is smoother, 
exhibits a much smaller overall volume distortion and shear (as visualized by 
the difference in principal stretches), and sharpens the compression-dilation 
dipoles around each single dislocation, as shown in Figure[6j 

The noisy PFC image of Figure [5] (right) is generated by adding 50% 
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Figure 6: Results of the image in Figure[5](left). From panel (a) to (d): crystal orientation, 
difference in principal stretches, curl of the inverse deformation gradient, and volume 
distortion. In each panel, the left figure shows the initial results from SST and the right 
one shows the optimized results from the variational method. Particularly in (c)-right, 
the average curl on each connected defect region is shown. 



Figure 7: Results of Figure[5] (right) using the same visualization as in Figure[6] 


23 






Figure 8: A comparison of the results obtained with different downsampling rates. From 
left to right, the synchrosqueezed transform downsamples the original image by a factor 
of 2, 4, and 8, and the runtime is about 50 minutes, 13 minutes and 20 seconds. Top 
panel: initial volume distortion provided by the synchrosqueezed transform. Bottom panel: 
optimized volume distortion. 
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Figure 9: Left: a TEM-image of GaN (courtesy of David M. Tricker, Department of 
Material Science and Metallurgy at the University of Cambridge). Right: a photograph 
of a bubble raft (courtesy of Barrie S. H. Royce at Princeton University). 


Gaussian random noise. Obviously, this leads to strong artifacts in the es¬ 
timated deformation Gq. Remarkably, after the optimization we retrieve an 
estimate G almost as good as in the noiseless case as shown in Figure[7| which 
demonstrates the robustness of our method. 

To show that the performance of onr method is not really depending 
on the mesh size chosen for the optimization algorithm, we compare the 
results obtained with different downsampling rates in the synchrosqueezed 
transform. Figure [8] summarizes the comparison for the PFC example. Even 
though the grain boundary geometry is complicated in this example, our 
algorithm obtains similar results for different downsampling rates. 

3.2. Experimental crystal images 

The first real crystal image, a TEM-image of GaN of size 420 x 444 
(Figure [9] left), contains a string of dislocations forming a large angle grain 
boundary. The artificial strong spatial variation of crystal orientation, shear 
and volume distortion in the SST result is greatly reduced after applying the 
optimization as shown in FigureflO} Even more importantly, the unphysical 
curl away from the defects is completely removed such that the curl of G is 
fully concentrated in the single defect regions around each dislocation (the 
total curl of each region equalling the dislocation’s Burgers vector). 


25 











Figure 10: Results of crystal analysis for Figure[9] (left), using the same visualization as in 
Figurc[6] 



Figure 11: Results of crystal analysis for Figure[9] (right), using the same visualization as 
in Figurc[6j 
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Figure 12: A TEM-image of A1 (courtesy of the National Center for Electron Microscopy 
at the Lawrence Berkeley National Laboratory). 


The second example is a photograph of a bubble raft with strongly dis- 
ordered and blurry grain boundaries (Figure |9| r ight). The image size is 
223 x 223 pixels. The SST result (see Figure 11) shows a spurious strong 
shear of the local crystal structure close to the grain boundaries, especially 
near the triple junction. One of the reasons for this behavior is that the SST, 
like any wavelet type transform, extracts directional and strain information 
from image patches, which here have to be larger than the unit cells. Thus, 
the grain boundaries are diffused, and information near the grain boundary 
is not trustworthy. The optimization can mostly remedies this effect. 

The last experimental example is a TEM-image of a twin and a high 
angle grain boundary in A1 (Figure [l2|). The image size is 424 x 634 pixels. 
The crystal structures in each grain are also slightly stretched. Although 
this example is very challenging, the SST-based analysis can still provide 
an accurate defect region estimate and a reasonable initial guess Go (see 
Figure[l3|). After optimization, we obtain a curl-free inverse deformation 
gradient G in the grain interior. The difference of principal stretches becomes 
smaller and the volume distortion gets closer to zero outside the defect region. 


4. Discussion 

We have presented a variational method to improve an initial analysis 
of a crystal image provided by the SST-based method developed in [ TtTj . 
By minimizing the elastic energy of the local inverse deformation gradient 
outside the SST-estimated defect region, we obtain an optimized deformation 
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Figure 13: Results of crystal analysis for Figure[l2] using the same visualization as in 
Figure[6] 

gradient G that better agrees with physical understanding. The desired 
information about the local crystal state is immediately available from G: 
the polar decomposition of G gives crystal orientations; the determinant of 
G indicates volume dilation and compression; the curl of G indicates defect 
regions, etc. 

The variational method in this paper presents a first step towards regu¬ 
larizing the local inverse deformation gradient provided by SST. In this step, 
we only focus on the optimization outside defect regions. In a next step, 
one has to investigate variational models to also obtain a meaningful esti¬ 
mate of the deformation gradient inside defect regions, which would provide 
a more accurate localization of defects. One possible solution could be a 
nested optimization model with an inner and an outer loop, where the inner 
loop performs the optimization outside the defect regions as proposed in this 
paper. The outer loop then tries to improve curl G inside the defect regions. 
We will leave this to future works. 
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